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HIGHLIGHTS 


►  Two  computational  models  of  the  air-cathode  of  a  PEMFC  are  developed  using  a  two-phase  approach. 

►  The  novelty  of  the  models  is  the  application  of  experimentally  determined  properties  of  the  GDL. 

►  Common  empirical  correlations  for  GDL  behavior  are  shown  to  underpredict  saturation. 

►  Experiments  show  that  the  model  can  predict  performance  based  on  measured  GDL  characteristics. 
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In  this  work  two  models  of  the  air-cathode  of  a  PEM  fuel  cell  are  proposed  using  a  two-phase  non¬ 
equilibrium  approach.  The  novelty  of  both  models  is  that  they  use  experimentally  determined  porosity, 
capillary  pressure  relationships,  and  permeability  of  the  gas  diffusion  material  in  order  to  demonstrate 
whether  or  not  the  use  of  these  properties  enhances  the  accuracy  of  PEM  fuel  cell  models  to  predict 
performance  curves  and  liquid  water  saturation  distribution.  The  first  model  is  a  1-D  model,  where  the 
effect  of  using  experimental  water  transport  properties  of  the  GDL  is  assessed  in  terms  of  water  satu¬ 
ration  distribution  in  the  GDL.  The  second  model  is  a  2-D  model  used  to  predict  experimental  polari¬ 
zation  curves.  It  was  found  that  the  implementation  of  the  experimentally  determined  water  transport 
characteristics  of  the  GDL  predicts  more  liquid  water  saturation  than  using  empirical  correlations  for  the 
capillary  pressure  curves  and  permeability.  It  was  also  observed  that  polarization  curves  of  two  cells 
using  different  GDL  material  can  be  predicted  accurately  using  the  appropriate  GDL  properties. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

PEM  fuel  cell  modeling  has  been  employed  for  many  years,  as  an 
effective  way  to  design,  improve,  and  predict  the  performance  of 
these  devices.  Accurate  modeling  is  not  an  easy  task  due  to  the  two- 
phase  phenomena  exhibited  by  operating  PEM  fuel  cells,  namely 
the  liquid  water  transport  and  the  reactant  gas  transport  through 
the  porous  GDL  There  are  three  approaches  to  study  the  two-phase 
flow  in  the  porous  media  of  the  fuel  cell:  (1)  the  pore  network 
approach  is  based  on  stochastic  procedures  to  generate  a  network 
of  pores,  aiming  to  reproduce  the  intricate  matrix  of  a  GDL  to 
predict  the  functional  relationship  of  the  water  transport  charac¬ 
teristics  [1-4];  (2)  the  direct  method  which  is  based  on  a  micro¬ 
scopic  treatment  of  the  GDL  in  which  the  actual  structure  of 
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the  porous  media  is  reconstructed  by  microscopic  imaging  [5];  and 
finally,  (3)  the  continuum  medium  approach  which  is  widely 
applied  in  the  literature  and  has  many  advantages  over  the  other 
methodologies,  including  the  capability  to  provide  a  complete  cell- 
level  model  and  low  computing  time  compared  to  higher  fidelity 
approaches. 

Continuum  models  refer  to  the  traditional  models  broadly  used 
in  many  applications  involving  a  flow  through  a  porous  media; 
some  applications  include  petroleum  engineering  and  subsurface 
hydrology.  The  porous  medium  is  treated  as  a  hypothetical  effective 
continuum  and  the  models  involve  volume-averaged  quantities 
such  as  saturation  and  rely  on  phenomenological  relationships 
such  as  the  generalized  Darcy’s  law.  Key  parameters  in  this 
approach  are  the  permeability,  relative  permeability  and  the 
capillary  pressure-saturation  curve  [1]. 

Two  main  modeling  schemes  can  be  distinguished  in  the 
continuum  models.  The  first  one  was  proposed  by  Wang  and 
Chen  [6].  These  authors  presented  a  novel  formulation  of  the 
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transport  processes  based  on  a  multicomponent  system  for  porous 
media.  The  model  equations  developed  include  mass  conservation, 
momentum  conservation,  species  transport,  and  the  energy  equa¬ 
tion.  The  multi-phase  mixture  model  is  commonly  referred  as  the  M2 
model.  One  of  the  main  features  of  the  M2  model  is  that  the  two- 
phase  stream  (gas-liquid)  is  treated  as  a  single  stream  by  introduc¬ 
ing  averaged  mixture  properties  and  variables  in  the  governing 
equations.  This  model  has  been  widely  employed,  from  its  initial 
development  [7]  on  transport  of  organic  compounds,  until  more 
recent  applications  to  PEM  fuel  cells  two-phase  models  [8-17]. 

Alternatively  to  the  M2  method,  there  is  a  simpler  model  often 
referred  as  the  multi-fluid-model  or  MFM.  This  method  is  based  on 
solving  the  conservation  equations  for  the  liquid  phase  and  the  gas 
phase  separately  using  source  and  sink  terms  to  account  for  mass 
exchange  among  phases.  The  MFM  has  been  widely  reported  in  the 
literature  when  dealing  with  two-phase  flow  analysis  in  PEMFCs 
[18—35].  In  the  MFM,  Darcy’s  equation  is  applied  in  the  GDL  to  each 
phase.  It  is  usually  assumed  that  the  gas  phase  pressure  gradient  is 
negligible;  this  assumption  is  called  the  unsaturated  flow  theory  or 
UFT.  The  UFT  has  been  applied  in  most  of  the  work  that  uses  the 
MFM  approach.  With  this  assumption,  the  gradient  of  the  capillary 
pressure  is  equal  to  the  negative  of  the  liquid  phase  pressure 
gradient: 


VPC  =  VPg  -  VPj  =  -VPj  (1) 

as  a  result,  the  liquid  phase  transport  equation  in  the  GDL  is 
transformed  to: 

v  (ftt)  =  v  (-'Kvp])  -v.(ffivpc)  .0  (2) 

where  k  is  the  effective  permeability  of  the  liquid  phase,  namely, 
the  absolute  permeability,  I<  (a  property  of  the  porous  media),  times 
the  relative  permeability  of  the  liquid  phase,  kT \  (a  function  of  the 
flow  conditions,  the  porous  matrix,  the  saturation  of  the  pores  with 
liquid  water,  etc.).  Using  the  chain  rule,  the  gradient  of  the  capillary 
pressure  can  be  expressed  as  a  function  of  saturation  knowing  that 
there  is  a  functional  relationship  Pc(s): 
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where  the  quantity  inside  the  inner  parenthesis  is  often  called  the 
liquid  water  diffusion  coefficient  Dw. 

So  far,  the  two  common  approaches  applied  to  the  analysis  of 
two-phase  flow  in  PEMFCs  have  been  briefly  described.  A 
comparison  between  these  two  approaches  was  addressed  by 


Pasaogullari  and  Wang  [9].  They  presented  a  ID  model  of  the 
cathode  side  GDL  of  a  PEMFC  where  the  flow  channels  and  the 
catalyst  layer  were  considered  as  boundary  conditions.  The  only 
equation  solved  in  [9]  was  the  liquid  mass  transport  equation 
using  both  approaches,  UFT  and  M2  model.  Pasaougullari  and 
Wang  found  that  the  UFT  under  predicts  saturation  in  comparison 
with  the  M2  approach;  they,  also  found  that  the  oxygen  distri¬ 
bution  in  the  GDL  was  also  underpredicted  by  the  UFT.  The  water 
saturation  obtained  by  the  UFT  was  slightly  below  the  M2 
prediction;  on  the  other  hand,  the  oxygen  distribution  did  show 
a  considerable  difference.  The  pressure  distribution  of  the  gas 
phase  was  compared  using  both  methods  and  a  deviation  of  2  kPa 
at  2  A  cm-2  in  a  GDL  of  0.3  mm  thick  was  observed. 

Regardless  of  the  approach  used  to  simulate  water  transport 
through  the  porous  media  of  a  PEM  fuel  cell,  two  constitutive 
relationships  are  necessary:  the  capillary  pressure-saturation  rela¬ 
tionship  and  the  liquid  phase  relative  permeability-saturation 
relationship.  These  expressions  are  necessary  due  to  the  capillary 
transport  mechanism  of  the  liquid  water  in  the  GDL.  Table  1  shows 
some  of  the  first  assumptions  found  in  the  literature  for  the  func¬ 
tional  relationships  required  for  a  continuum  model.  He  et  al.  [18] 
just  assumed  a  constant  liquid  water  diffusion  coefficient  and 
a  linear  function  for  the  relative  permeability  of  the  liquid  water. 
Nartarajan  and  Nguyen  [19]  proposed  an  expression  for  the 
gradient  of  the  capillary  pressure  curve  where  the  constants  were 
adjusted  via  experimental  polarization  curves  and  the 
relative  permeability  was  defined  as  a  linear  function  of  saturation 
Wang  et  al.  [8]  were  one  of  the  first  groups  to  use  a  modification 
of  the  empirical  Leverret  function  for  well-sorted  sand  to  evaluate 
the  capillary  pressure-saturation  relationship  in  a  PEM  fuel  cell 
GDL.  These  empirical  relationships  are  written  in  terms  of  reduced 
saturation  (S)  instead  of  saturation  (s).  In  order  to  simplify  the 
modeling  efforts,  Nam  and  Kaviany  evaluated  the  slope  of  the 
Leverret  function  in  the  range  0  <  S  <  0.5  to  represent  the  capillary 
pressure  curve  relationship.  Pasaogullari  et  al.  [12],  proposed  a 
modified  form  of  the  Leverret  function  to  describe  the  effects  of 
having  a  hydrophobic  or  hydrophilic  porous  media. 

The  modification  of  the  empirical  Leverret  correlation  for  well- 
sorted  sand  proposed  by  Pasaogullari  et  al.  [12]  and  the  liquid 
phase  relative  permeability  empirical  relationship  for  the  same 
material  have  been  widely  used  among  researchers  involved  in  the 
modeling  of  PEM  fuel  cells.  The  popularity  of  this  method  can  be 
seen  in  the  large  number  of  modeling  articles  referring  to  the 
Leverret  function  and  the  cubic  expression  for  relative  perme¬ 
ability  of  the  liquid  phase  as  the  most  appropriate  or  simplest 
manner  to  address  the  effects  of  capillary  motion  of  liquid  water 
through  the  GDL.  The  disadvantage  to  this  approach  is  that  it  is  not 
well-connected  to  measurable  characteristics  of  the  GDL  material 


Table  1 

Different  expressions  found  in  the  literature  for  capillary  pressure  and  liquid  water  relative  permeability  relationships. 


Pc(s)  relationship 

/<ri  (liquid  phase 
relative  permeability) 

Reference 

Dw  =  l  x  10-4  cm2  s”1 

kTi  =  s 

[18] 

krl  =  s  +  0.01 

[19] 

acosd 417  2.12(1  S)2  +  1.236(1  S)3] 

(I</e)V2 

fen  =  S3 

[8] 

d  Pr 

=  30321  [Pa]  Slope  of  Leverret  function  0  <  S  <  0.5 

fen  =  S3 

[20] 

pc=  ,C0S,!2i(s) 

(I</e)V2 

(  1.417S-2.12S2  +  1.262S3  9  >  90° 

J(s)  =  <  ,  o" 

1 1.417(1  -S)-  2.12(1  -  S)2+l. 262(1  —  S)3  9  <  90° 

fen  =  S3 

[12] 
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and  hence  cannot  provide  insight  into  how  GDL  characteristics 
affect  performance.  In  this  work,  the  experimental  characteriza¬ 
tion  performed  on  carbon  paper  Toray  090  by  Ramos-Alvarado 
et  al.  [36]  is  taken  as  the  base  of  two  numerical  models  of  the 
air-cathode  of  a  PEM  fuel  cell.  First,  a  1-D  model  is  developed 
aiming  to  show  that  the  empirical  Leverret  function  and  relative 
permeability  correlation  for  well-sorted  sand  are  inappropriate 
when  liquid  water  saturation  distribution  is  predicted  for  carbon 
paper  under  different  wet-proof  treatments.  The  second  model  is 
a  2-D  representation  of  the  cathode  side  of  a  PEM  fuel  cell  oper¬ 
ating  at  very  high  stoichiometric  conditions.  The  objective  of  the 
2-D  model  is  to  predict  polarization  using  different  GDL  materials 
with  different  water  transport  properties  and  compares  the 
predictions  with  experimental  polarization  curves. 

2.  Model  development 

In  this  section,  the  modeling  efforts  followed  to  assess  the 
impact  of  using  experimentally  determined  water  transport  char¬ 
acteristics  to  predict  liquid  water  saturation  and  polarization  in 
a  PEMFC  are  presented.  First,  a  1-D  model  of  the  air-cathode  of 
a  PEM  fuel  cell  is  developed  to  study  the  distribution  of  liquid  water 
saturation.  This  model  was  run  using  the  widely  used  empirical 
correlations  for  well-sorted  sand  and  with  the  experimentally 
determined  capillary  pressure  curves  and  relative  permeability 
functions  obtained  in  [36].  Secondly,  a  2-D  model  of  the  air-cathode 
of  a  PEMFC  is  developed  to  predict  experimental  polarization 
behavior.  This  2-D  model  is  run  using  experimentally  determined 
water  transport  characteristics  and  electrochemical  parameters 
including  exchange  current  density,  open  circuit  voltage,  and  the 
transfer  coefficient. 

2.1.  Modeling  assumptions 

The  simplifications  implemented  to  the  1-D  model  are  listed 
below. 

1.  Steady  state  operation  of  the  fuel  cell. 

2.  Isothermal  operation. 

3.  Liquid  water  is  introduced  at  the  catalyst  layer-GDL  interface 
due  to  the  cathode  electrochemical  reaction  and  due  to  electro- 
osmotic  drag  from  the  anode. 

4.  The  gas  phase  pressure  is  constant,  which  means  a  null  gas 
phase  pressure  gradient  throughout  the  domain.  The  UFT 
assumption  is  applied. 

5.  No  irreducible  saturation  is  present  in  the  gas  diffusion  media. 

6.  There  is  no  equilibrium  of  phases  between  liquid  water  and 
water  vapor  (i.e.  water  can  move  between  phases). 

7.  A  1 D  model  approach  is  appropriate  to  compare  the  effects  of 
using  different  constitutive  relations  in  the  prediction  of  liquid 
water  saturation  distribution. 

For  the  development  of  the  2-D  model,  the  assumptions  from  1 
to  6  were  taken,  except  that  the  2-D  model  considers  a  non- 
isothermal  operation.  The  2-D  analysis  is  based  on  reactant  and 
product  concentrations  at  a  particular  point  along  the  flow  channel 
and  assumes  that  the  change  in  composition  along  the  flow  channel 
is  negligible  (i.e.  operation  with  a  high  stoichiometric  ratio). 

2.2.  Computational  domains 

Fig.  1  shows  the  cross-section  of  a  PEMFC.  Both  sides  compo¬ 
nents,  anode  and  cathode,  are  shown.  The  effect  of  compression  of 
the  current  collector  over  the  gas  diffusion  media  is  depicted,  as 
well  as  the  catalyst  layers  and  the  polymeric  electrolyte  membrane. 


Fig.  1.  Cross-section  of  a  PEMFC 


The  models  developed  will  study  only  the  cathode  side  of  the  fuel 
cell  due  to  the  higher  activation  barrier  of  the  cathode  reactions  in 
comparison  with  the  anode  reactions.  The  computational  domains 
will  be  described  in  the  following  paragraphs. 

The  computational  domain  of  the  1-D  model  consists  of  the  one¬ 
dimensional  length  through  plane  of  the  cathode  GDL.  Fig.  2  shows 
the  geometry  studied.  The  boundaries  of  this  model  are  the  inter¬ 
face  with  the  catalyst  layer  (x  =  0)  and  the  interface  with  the  gas 
flow  channel  (x  =  L).  In  this  model  L  =  287  pm,  which  is  the  average 
thickness  of  the  gas  diffusion  media  studied. 

In  order  to  simplify  the  numerical  computations  in  the  2-D 
model,  a  discrete  region  is  taken  as  computational  domain;  see 
Fig.  3.  The  computational  domain  consists  of  two  subdomains:  (1) 
Di  is  the  PEM  and  (2)  D2  is  the  GDL,  both  subdomains  are  half  of 
the  cross-section  shown  in  Fig.  1.  In  this  computational  domain, 
the  catalyst  layer  was  reduced  to  an  interface  between  sub- 
domains  D\  and  D2.  The  effects  of  the  catalyst  layer  will  be 
considered  in  the  electrochemistry  of  this  model.  The  dimensions 
shown  in  Fig.  3  are:  (1)  the  uncompressed  thickness  of  the  gas 
diffusion  layer  under  the  gas  flow  channels  (xu);  (2)  the 
compressed  thickness  of  the  gas  diffusion  layer  under  the  collector 
plate  (xc);  (3)  the  in  plane  length  of  the  uncompressed  GDL  (yu); 
(4)  the  in  plane  length  of  the  compressed  portion  of  the  GDL  (yc); 
and  (5)  the  PEM  thickness  (xm). 

2.3.  Governing  equations 

Table  2  shows  the  governing  considered  in  the  analysis  as  well 
as  information  regarding  the  application  of  such  equations. 

Given  that  air  is  the  reactant  gas  in  the  cathode,  Eq.  (4)  repre¬ 
sents  three  equations,  one  for  each  component:  02,  H20,  and  N2; 
however,  only  two  equations  are  solved,  02  and  H20,  and  the  mass 
fraction  of  N2  is  calculated  by: 

X>  =  1  O) 

1  =  1 

The  source  term  in  Eq.  (4)  accounts  for  species  consumption  or 
generation  in  the  computational  domain.  In  this  model,  only  water 
goes  through  phase  change  and  is  modeled  by  [18]: 


1— 

— > 

Catalyst 

X 

Gas  flow 

layer 

GDL 

Fig.  2.  1-D  model  computational  domain. 
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In  order  to  ensure  conservation  of  mass  between  the  liquid 
water  evaporating  (from  the  liquid  phase)  and  the  water  vapor 
gained  (to  the  gas  phase)  or  the  liquid  water  condensating  (to  the 
liquid  phase)  and  the  water  vapor  lost  (from  the  gas  phase),  the 
source  term  of  Eq.  (5)  is  defined  as: 


ri  =  ~ru2o  (13) 

The  diffusivity  of  species  was  calculated  using  a  symmetric 
multicomponent  Fick  model.  In  order  to  use  this  multicomponent 
model,  the  Maxwell— Stefan  binary  diffusivity  is  determined  by: 


kpP-75  ( 
(,V3  + ,V3  )2VMi 


(14) 


Fig.  3.  2-D  model  computational  domain,  dimensions  and  subdomains. 


rH20  =  /<c 


e(\  S)Xh2o  ( „  n  Dsat 


RT 

£Sp[ 


Mi 


h2o 


xh2oP-Ph2o)  Q 


Xh2oP-Ph2o)(1-<J) 


(10) 


where  /<d  is  a  constant  equal  to  3.16  x  10  8  Pa  m2/s,  T is  the  absolute 
temperature  in  K,  M,  is  the  molar  mass  of  species  i  in  kg/mol,  P  is  the 
pressure  expressed  in  Pa,  and  u;  equals  the  molar  diffusion  volume 
of  species  i  expressed  in  m3/mol.  The  effective  Maxwell— Stefan 
diffusivities  are  calculated  as  suggested  in  [20]: 


Df  =  D« 


e-0.11 
1  -0.11 


0.785-1 


(1-S) 


(15) 


where  kc  and  ke  are  the  condensation  and  evaporation  constants,  s 
is  the  local  liquid  water  saturation,  P  is  the  pressure  of  the  mixture, 
Pj^o  is  the  saturation  pressure,  T  is  the  absolute  temperature  in  K, 
xh20  is  the  molar  fraction  of  water  vapor,  Mh2o  is  the  water  molar 
mass,  pi  is  the  liquid  water  density,  and  q  is  the  switch  function  that 
allows  to  change  between  condensation  and  evaporation: 


=  1_  I^H2oP__^o^ 

2  +  2  xH2oP-P%0 


(11) 


Water  vapor  is  either  consumed  or  generated  due  to  phase 
change  in  the  air-cathode  of  a  PEMFC.  Oxygen  is  consumed  at  the 
catalyst  layer  interface  only.  Nitrogen,  on  the  other  hand,  does  not 
react  at  all  in  the  computational  domain;  thus,  a  flow  stagnation 
condition  can  be  established  in  this  component,  allowing  to 
calculate  the  bulk  velocity  of  the  gas  mixture  as: 


Table  2 

Governing  equations. 


Equation 

Mathematical  expression 

1-D 

Model 

2-D 

Model 

Species 

transport 

(  n  M  \ 

V-  1  pso>iH  -  pgu,  ^Dg-^ViOj  J 

=  n  (4) 

V* 

Liquid 

water 

(5) 

V* 

transport 

Energy 

V(VgCpuT  kfvT)  =  rT 

(6) 

Electric 

potential 

field 

V-(<7  eV0e)  =  0 

(7) 

V* 

Protonic 

potential 

field 

V-(<T;V0j)  =  0 

(8) 

(Di) 

In  Eq.  (15),  the  term  inside  the  square  brackets  accounts  for  the 
effect  of  the  reduction  of  diffusion  due  to  porosity  of  the  GDL,  and 
the  second  term  accounts  for  the  reduction  of  diffusivity  due  to 
a  partially  saturated  GDL  with  liquid  water.  At  the  end,  the  effective 
diffusivities  are  arranged  in  a  symmetric  matrix  used  to  calculate 
the  multicomponent  Fick  diffusivities. 

In  the  energy  equation,  the  source  term,  rT,  represents  the 
thermal  effect  of  phase  change,  and  is  calculated  as: 


rT  =  rH2oAftfg  (16) 

where  Ahfg  is  the  water  enthalpy  of  phase  change.  In  the  energy 
equation,  the  term  kjf  represents  the  effective  thermal  conduc¬ 
tivity  of  the  computational  domain,  namely,  carbon,  liquid  water, 
and  humid  air: 

kf  =  (1  -  e)kcarbon  +  eski  +  e(l  -  s)kg  (17) 

The  current  density  is  calculated  using: 

j  =  jo(l  -  s)  ^exp  4act)  (18) 

where  jo  is  the  exchange  current  density,  C02  is  the  concentration  of 
oxygen  at  the  catalyst  layer  interface,  cgj  is  the  reference  concen¬ 
tration  of  oxygen,  Pis  Faraday’s  constant,  a  is  the  transfer  coefficient 
and  relates  to  the  symmetry  of  the  reaction,  n  is  the  number  of 
electrons  transferred  in  the  electrochemical  reactions,  and  ^ac t  is 
the  activation  overpotential  given  by: 

Vact  =  Voc  —  ^cell  —  ^ohm  (19) 

where  Voc  is  the  reversible  cell  voltage  or  open  circuit  voltage,  VCeii 
is  the  actual  cell  voltage,  and  ^hm  is  the  ohmic  overpotential 
which  is  calculated  as  the  summation  of  the  GDL  and  membrane 
loses: 
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Vohm  =  2r?GDL  +  ^PEM  (20) 

where  the  GDL  overpotential  is  assumed  to  be  similar  for  anode 
and  cathode  side,  tjgdl  is  calculated  as  the  electric  potential  drop 
and  ?7pem  is  the  protonic  potential  drop  in  the  PEM.  It  can  be  seen 
that  ?7ohm  includes  the  effect  of  the  PEM  overpotential  only  in  the 
2-D  model. 

2.4.  Boundary  conditions 

The  boundary  conditions  applied  to  the  1-D  model  are 
summarized  in  Table  3. 

Fig.  4  shows  all  the  boundaries  defined  in  the  2-D  model 
developed  where  Bi  is  the  catalyst  layer  reduced  to  an  interface,  B2 
is  the  interface  with  the  gas  flow  channel,  B3  and  B4  are  the  inter¬ 
faces  with  the  current  collector  plate,  B5  is  the  symmetry  plane  at 
the  midpoint  of  the  collector  shoulder,  B6  is  the  symmetry  plane  at 
the  midpoint  of  the  flow  channel,  and  B7  the  interface  between  the 
anode  side  of  the  fuel  cell  and  the  PEM.  Table  4  summarizes  the 
boundary  conditions  set  at  each  boundary. 

Symmetry  conditions  where  defined  at  boundaries  B5  and  B6.  At 
B7  the  protonic  potential  was  defined  as  the  corresponding  anode 
potential,  namely  </>i  =  0.  At  B2,  the  mass  fractions  w02,gfc  and 
wh2o,gfc  are  calculated  based  on  air  relative  humidity  supplied  as 
reactant  to  the  cathode. 

2.5.  Parameters  and  operating  conditions  of  the  model 

The  parameters  employed  in  the  numerical  model  are  listed  in 
Table  5. 

The  properties  of  carbon  paper  Toray  090  were  taken  from 
Ref.  [36].  Two  capillary  pressure  curves,  one  for  non-wet-proofed 
material  (untreated)  and  one  for  a  GDL  material  with  20%  of  PTFE 
were  fitted  as  shown  in  Fig.  5. 

Table  6  presents  the  properties  of  the  GDL  material  that  was 
characterized  in  [36],  properties  that  are  going  to  be  used  in  the 
numerical  models. 

2.6.  Numerical  procedures 

COMSOL’s  default  solver  was  used  for  the  solution  of  the  gov¬ 
erning  equations.  A  direct  method  was  selected  to  solve  the 
stationary  problems  (1-D  and  2-D)  using  the  MUMPS  solver.  The 
parallel  sparse  direct  linear  solver  MUMPS  (MUltifrontal  Massively 
Parallel  Sparse  Direct  Solver)  works  on  general  systems  of  the  form 
Ax  =  b.  MUMPS  uses  several  preordering  algorithms  to  permute  the 
columns  and  thereby  minimizing  the  fill-in.  MUMPS  is  multi¬ 
threaded  on  platforms  that  support  multithreading  and  also 
supports  solving  on  distributed  memory  architectures  through  the 
use  of  MPI.  MUMPS  also  includes  out-of-core  capabilities.  The 
MUMPS  out-of-core  solver  stores  the  LU  factors  on  the  hard  drive. 


Table  3 

1-D  model  boundary  conditions 


Equation 

GDL— catalyst  layer  interface 

GDL— gas  flow  channels 
interface 

Oxygen  transport 

Water  vapor 
transport 

™o2  =  -Mo2Jp 

”1h2o  =  0 

1  -  wh2o 

W°2  “  1  ,  r 

1  +  rN2/02 

Mu  n  Pun 

^0  =  ™ 

Liquid  water 
transport 

Electric  potential 
field 

rh\  =  MHz  ojrpO  +  2anet) 

rnx  =  hm(s  —  sGfc) 

-n-((TeV0e)  =  -j 

0e  =  ^cell 

B 


b2 


Fig.  4.  Boundaries  of  the  2-D  model. 

This  minimizes  the  internal  memory  usage.  The  price  to  pay  is 
longer  solution  times  because  it  takes  longer  time  to  read  and  write 
to  disk  than  using  the  internal  memory  [37]. 

3.  Experimental  techniques 

A  commercial  Ion  Power  Inc.,  Nation  catalyst  coated  membrane 
was  used  as  MEA.  Both  sides  of  this  coated  membrane  have 
a  platinum  content  of  0.3  mg  cm-2.  Two  different  GDL  materials 
were  used:  Toray-090  untreated  (without  wet  proofing)  and  the 
same  material  with  a  PTFE  content  of  20%  (wet  proofing  20%  wt). 
The  gasket  material  used  was  PTFE  coated  fiber  glass.  A  single 
serpentine  pattern  served  as  gas  distributor.  The  total  length  of  the 
serpentine  is  322  mm  with  square-shaped  channels  of  1  mm.  The 
total  active  area  of  this  assembly  is  4.84  cm2. 

The  test  conditions  were  determined  so  that  a  high  stoichio¬ 
metric  flow  was  attained.  By  using  a  high  stoichiometric  ratio,  the 
depletion  effects  in  the  channels  are  not  that  important;  hence,  the 
2-D  mathematical  model  developed  is  suitable  for  describing  the 
cell  phenomena.  The  experimental  test  conditions  are  listed  in 
Table  7. 

Following  the  construction  of  the  test  cells,  these  were  installed 
in  a  commercial  fuel  cell  test  station  manufactured  by  Fuel  Cell 
Technologies  Inc. 

3.1.  Repeatability  of  experiments 

After  reaching  a  steady  state  operation  of  the  test  cells  under  the 
conditions  described  in  the  last  section,  the  polarization  curves 
were  obtained.  Three  different  experiments  were  conducted  by 
varying  the  recording  time  of  current  at  constant  voltage  condi¬ 
tions.  Recording  times  of  20,  30,  and  40  s  were  tested  and  are 
named  as  Experiment  1,  Experiment  2,  and  Experiment  3  respec¬ 
tively.  Fig.  6  shows  the  repeatability  achieved  during  the  experi¬ 
mentation  with  the  two  different  single  cells  assembled.  The 
polarization  curves  for  the  three  experiments  were  found  to  be  very 
similar,  in  each  cell.  Such  results  demonstrate  a  good  repeatability 
of  the  experimental  procedures. 

3.2.  Polarization  curves  and  Tafel  plot 

Fig.  7  shows  that  the  fuel  cell  using  an  untreated  (non-wet- 
proofed)  GDL  material  presents  a  higher  performance  than  the  fuel 


Table  4 

2-D  model  boundary  conditions. 
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Equation 

Bi 

b2 

B3,  B4 

Oxygen  transport 

™o2  =  -Mo2-L 

™02  =  ^d(wo2  -wo2,gfc) 

3- 

P 

II 

o 

Water  vapor  transport 

mH2o  =  0 

™H20  =  ^d(wH20  -  wH20,GFc) 

mH2o  =  0 

Liquid  water  transport 

mi  =  Mh2  ojrp(l  +  2dnet) 

m\  =  hm(s-sG¥C) 

rii\  =  0 

Energy 

-n-(pgcpuT  -  kfvi)  =  /T(Voc  -  Vce\\)j 

q"  =  hj(T  —  T0) 

H 

II 

Electric  potential  field 

— n-  (o"eV0e)  =j 

-n-(aeV<l>e)  =  0 

</>e  —  Vcen 

Protonic  potential  field 

-n-(^/V0i)  =  -j 

N/A 

N/A 

Table  5 

Numerical  model  parameters. 


Description 

Symbol 

Value 

Units 

Mass  ratio  of  N2  to  02 

rN2/02 

3.2917 

- 

Isobaric  gas  pressure 

Pg 

101.325 

kPa 

Cell  temperature 

Lee  11 

353 

K 

GDL  porosity 

£ 

Varies  with 
material 

Net  water  transfer  coefficients 

®net 

0.5 

- 

Condensation  constant 

kc 

100 

s-1 

Evaporation  constant 

ke 

1 

atm-1  s-1 

Mass  convection  coefficient 

hm 

1000 

kg  rrr2  s-1 

Heat  transfer  convection  coefficient 

hT 

1.13 

W  m2  K1 

Diffusive  convection  coefficient 

ho 

Varies  with 
species 

kg  rrr2  s-1 

Carbon  thermal  conductivity 

^carbon 

6.75 

W  m-1  K-1 

Liquid  water  thermal  conductivity 

ki 

0.58 

W  m-1  K-1 

Gas  phase  thermal  conductivity 

kg 

0.02 

W  m-1  K-1 

Specific  capacity  of  the  gas  phase 

cp 

1008 

J  kg"1 K1 

Heat  partition  factor 

/t 

0.7 

- 

cell  using  wet-proofed  material.  These  results  may  look  contro¬ 
versial  in  regards  to  what  is  expected  from  the  PTFE  treatment, 
namely,  a  better  water  management  in  the  cathode  side;  never¬ 
theless,  the  following  must  be  considered:  (1 )  it  was  demonstrated 
that  the  material  with  20%  of  PTFE  has  a  very  low  permeability  in 
comparison  with  a  compressed  non  wet-proof  material  [36];  such 
a  characteristic  leads  to  a  reduced  diffusion  of  fuel  in  the  anode 
side,  likewise,  it  leads  to  a  low  diffusion  of  oxygen  in  the  cathode 
side;  (2)  the  operating  conditions  used  in  the  experiments 
encourage  good  water  removal,  namely,  the  cathode  relative 


Fig.  5.  Curve  fitting  to  experimental  capillary  pressure  curves  taken  from  Ref.  [36]. 


Table  6 

Experimentally  determined  properties  of  carbon  paper  Toray  090  [36]. 

Material  Pc(s)  relationship  [Pa]  kT\  I<  [m2]  e 

dP 

Untreated  =  (2.7322  -  11.6345  +  11.487S2)  x  104  S5  8.236xl0~12  0.78 
dS 

20%  PTFE  ^  =  (5.3907- 21. 641S  +  20.8519S2)  x  104  S2  1.098xl0~12  0.73 


humidity  inlet  condition  was  70%;  (3)  the  high  frequency  resistance 
of  the  cells  was  measured  as  0.0968  Q  cm2  and  0.1113  Q  cm2,  for  the 
cell  using  non-wet-proofed  and  for  the  cell  using  wet-proofed  GDL 
respectively;  and  (4)  the  hydrophobic  characteristics  of  the  GDL 
reduces  the  amount  of  moisture  in  the  MEA  and  thus  its  conduc¬ 
tivity  is  reduced.  These  arguments  aid  to  explain  the  behavior 
observed  in  Fig.  7;  further,  there  are  references  in  the  technical 
literature  that  support  these  findings  [38-41]. 

Wet-proofing  the  GDLs  leads  to  a  better  water  removal  but  also 
decreases  the  gases  diffusivity  and  the  electric  conductivity  of  the 
material;  thus,  an  optimal  PTFE  content  is  expected  to  exist  to 
account  for  this  trade-off.  Flowever,  there  are  many  answers  to  this 
problem  in  the  literature.  Velayutham  [38]  found  an  optimal  value 
of  35%  PTFE  content  in  terms  of  peak  power  density;  however,  this 
experimental  work  was  conducted  only  in  the  micro-porous  layer. 
Park  et  al.  [39]  found  an  optimal  value  of  PTFE  of  10%  also  in  the 
micro-porous  layer.  Jie-Cheng  and  Chien-Kung  [40]  conducted 
a  wide  study  at  different  humidity  temperatures  and  at  different 
PTFE  loadings  in  the  GDL  and  in  the  MPL.  They  found  an  optimal 
content  of  40%  PTFE  in  the  GDL  and  30%  PTFE  in  the  MPL  at  low 
humidity  conditions.  Park  et  al.  [41  ]  conducted  a  study  addressing 
the  effect  of  PTFE  content  in  the  GDL  (Toray  090)  on  cell  perfor¬ 
mance  at  different  cell  and  gases  temperatures.  The  work  presented 
in  [41  ]  includes  operating  conditions  very  similar  to  those  used  in 
this  work.  Park  et  al.  determined  that,  at  operating  conditions  like 
the  ones  used  in  this  study,  the  optimal  PTFE  content  is  10%  fol¬ 
lowed  by  0%,  and  for  PTFE  contents  of  20%  and  up,  the  cell  perfor¬ 
mance  decreases.  This  short  literature  review  on  the  effect  of  wet¬ 
proofing  in  GDL  material  on  cell  performance  is  presented  with  the 
objective  to  show:  (1)  that  no  clear  trends  are  observed  in  regards 
to  the  effect  of  GDL  wet-proofing  and  cell  performance,  and  (2)  that 
there  are  authors  reporting  similar  trends  of  cell  performance  like 
the  ones  presented  here,  under  similar  operating  conditions  and 
using  the  same  GDL  material. 


Table  7 

Experimental  test  conditions. 


Parameter 

Value 

Units 

Anode  flow  rate 

225 

seem 

Cathode  flow  rate 

550 

seem 

Cell  temperature 

353 

K 

Anode  gas  temperature 

323 

K 

Cathode  gas  temperature 

343 

K 

Anode/cathode  back-pressure 

0 

Pa 

382 
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Repeatability  of  experiments,  polarization  curves,  (a)  Single  cell  with  an  untreated  GDL.  (b)  Single  cell  with  a  GDL  with  20%  of  PTFE. 


The  Tafel  plot  was  obtained  for  both  test  cells  by  calculating  the 
logarithm  of  low  current  densities  in  the  activation  region  of  the 
polarization  curves;  however,  the  plots  were  similar  due  to  the  fact 
that  similar  MEA  material  (from  the  same  stock)  was  used  in  both 
cells.  The  resulting  Tafel  plot  is  shown  in  Fig.  8.  The  electrochemical 
parameters  calculated  from  this  plot  are  j0  =  49.5Am-2  and 
an  =  1.223. 


4.  Results 


In  this  section,  the  1-D  modeling  efforts  are  presented.  This 
model  was  developed  aiming  to  assess  the  impact  of  using  exper¬ 
imentally  determined  capillary  pressure  functions  and  relative 
permeability  functions  on  water  saturation  profiles  across  the  gas 
diffusion  media.  In  order  to  conduct  this  comparison,  the  widely 
reported  empirical  Leverett  correlation  and  the  empirical  correla¬ 
tion  for  liquid  water  relative  permeability  for  well-sorted  sand 
(/<ri  =  S3)  was  implemented  and  compared  to  the  model  using 
experimental  data. 

The  2-D  modeling  results  are  presented  in  this  section.  First, 
the  properties  of  a  compressed  non-wet-proofed  GDL  will  be  used 
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Fig.  7.  Polarization  curves  of  single  cells  using  different  GDL  material. 


in  the  numerical  model  to  see  if  polarization  can  be  predicted. 
Later,  the  properties  of  a  compressed  wet-proofed  GDL  are  going 
to  be  used  as  inputs  of  the  numerical  model  in  order  to  determine 
if  by  changing  these  properties,  polarization  can  be  predicted 
accurately. 

4.1.  1-D  model  results 

The  properties  of  a  compressed  non-wet-proofed  (untreated) 
material  were  implemented  into  a  1-D  model  using  the  properties 
and  operating  conditions  shown  in  Table  5  and  the  GDL  properties 
shown  in  Table  6.  A  different  1-D  model  was  studied  by  only 
changing  the  capillary  pressure  function  to  the  Leverett  function, 
see  Table  1  [12],  and  using  kr \  =  S3  to  describe  relative  permeability. 
The  saturation  profiles  obtained  are  compared  in  Fig.  9(a)  at  an 
operating  condition  of  1  A  cm-2. 

It  can  be  clearly  seen  that  the  empirical  correlations  of  the  type 
that  have  been  widely  used  in  the  literature  underpredict  satura¬ 
tion  in  the  GDL.  In  general,  by  using  the  experimental  water 
transport  characteristics,  the  water  saturation  prediction  is  nearly 


Fig.  8.  Tafel  plot. 
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Fig.  9.  Prediction  of  saturation  profiles  in  the  GDL.  (a)  Comparison  between  experimental  and  empirical  correlations,  (b)  Comparison  using  experimental  correlations  for  untreated 
and  wet-proofed  GDL 


twice  than  the  water  saturation  predicted  using  empirical  corre¬ 
lations.  It  can  be  observed  that  the  water  saturation  is  greater  near 
the  catalyst  layer.  This  is  due  to  the  liquid  water  generation  by  the 
electrochemical  reactions.  The  saturation  profiles  show  a  decrease 
near  the  gas  flow  channels  where  water  is  removed  by  convective 
mechanisms. 

The  effect  of  wet-proofing  on  the  saturation  profiles  is  depicted 
in  Fig.  9(b).  The  effect  of  a  wet-proofing  treatment  with  of  20% 
(PTFE)  is  remarkable.  Fig.  9(b)  shows  that  the  maximum  saturation 
using  the  material  with  20%  of  PTFE  exhibits  a  maximum  saturation 
of  2.4%  whereas  an  untreated  material  shows  a  maximum  satura¬ 
tion  of  14.6%.  Such  an  effect  is  attributed  to  the  hydrophobic 
characteristics  given  to  the  GDL  material  by  the  PTFE  treatment 
(wet-proofing). 

A  parametric  study  of  liquid  water  saturation  profiles  on 
untreated  compressed  GDL  is  presented  in  Fig.  10.  This  study  is 
conducted  in  order  to  assess  the  reliability  of  the  response  of  the 
model  to  different  operating  conditions.  Fig.  10(a)  is  a  parametric 
study  of  the  effect  of  the  inlet  air  relative  humidity  on  water 


saturation  profiles  at  an  operating  condition  of  1.4  A  cm-2.  It  is 
observed  that  as  the  relative  humidity  condition  in  the  cathode 
increases,  the  saturation  increases  due  to  the  high  water  content 
in  the  air  stream.  This  high  water  concentration  in  the  air  leads 
to  a  reduced  evaporation  rate  and  condensation  is  more  likely  to 
occur.  On  the  other  hand,  as  shown  later  in  this  section,  low 
humidity  conditions  drive  a  larger  evaporation  rate.  Fig.  10(b) 
shows  a  parametric  study  of  water  saturation  distribution  at 
a  constant  relative  humidity  condition  of  95%  and  a  ranges  of 
current  density  values.  It  can  be  observed  in  Fig.  10(b)  that  at 
larger  current  densities,  the  water  saturation  level  increases.  Such 
a  behavior  is  obvious  because  the  water  production  is  directly 
proportional  to  the  current  density.  By  comparing  Fig.  10(a)  and 
(b)  it  can  be  inferred  that  the  operating  current  density  operating 
condition  has  a  greater  impact  on  water  saturation  than  that  of 
the  inlet  air  relative  humidity. 

Fig.  11  presents  a  parametric  study  of  oxygen  distribution  in 
a  compressed  and  untreated  GDL.  Fig.  11(a)  shows  a  study  of  the 
effect  of  inlet  air  relative  humidity  on  oxygen  distribution  at 


Fig.  10.  Parametric  study  of  liquid  water  saturation  profiles  in  the  GDL  (a)  Cell  operating  at  1.4  A  cm  2.  (b)  Air  inlet  relative  humidity  of  95%. 
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Fig.  11.  Parametric  study  of  oxygen  distribution  in  the  GDL.  (a)  Cell  operating  at  1.4  A  cm  2.  (b)  Air  inlet  relative  humidity  of  95%. 


a  constant  current  density.  Relative  humidity  plays  an  important 
role  when  determining  the  concentration  of  oxygen  and  water 
vapor  in  humid  air.  Having  a  high  relative  humidity  condition 
means  a  lower  concentration  of  oxygen  in  humid  air;  conversely, 
a  low  humidity  represents  mixture  with  a  higher  concentration  of 
oxygen.  Fig.  11(b)  shows  a  parametric  study  of  oxygen  concentra¬ 
tion  at  a  constant  inlet  air  relative  humidity  condition  of  95%  and 
varying  current  density.  As  current  density  increases,  the 
consumption  of  oxygen  increases  and  thus  the  concentration 
gradient  which  drives  diffusion  must  increase.  It  can  be  seen  that  at 
low  current  densities  like  0.16  A  cnrr2  the  gradient  of  oxygen  across 
the  GDL  is  not  that  considerable  as  that  of  an  operating  condition  of 
2.1  A  cm-2. 

Whereas  oxygen  is  consumed  at  the  interface  between  the 
catalyst  layer  and  the  GDL,  this  interface  is  impermeable  to  water 
vapor.  Fig.  12(a)  shows  the  effect  of  inlet  air  relative  humidity  on 
water  vapor  distribution.  All  the  curves  shown  have  a  different 
starting  point  at  the  gas  flow  channels  interface  that  corresponds  to 
the  inlet  humidity  condition  defined  for  each  case.  It  is  observed 
that  water  vapor  mass  fraction  increases  in  direction  from  the  gas 
flow  channels  interface  towards  the  interface  with  the  catalyst 


layer.  The  increase  in  water  vapor  mass  fraction  is  due  to  the 
evaporation  that  takes  place  inside  the  GDL  as  predicted  by  the 
phase  change  model  used  (see  Eq.  (10)).  The  driving  mechanism  of 
phase  change  in  this  model  is  the  non-equilibrium  between  phases. 
A  measure  of  this  non-equilibrium  is  given  by  the  difference 
between  the  local  water  vapor  pressure  and  the  saturated  vapor 
pressure  at  the  local  temperature.  It  can  be  seen  how  the  gain  in 
mass  fraction  is  larger  for  low  humidity  (and  hence  low  vapor 
pressure)  conditions. 

Fig.  12(b)  shows  the  effect  of  different  current  density  operating 
conditions  on  water  vapor  mass  fraction  distribution.  All  the  cases 
studied  share  a  common  value  at  the  interface  with  the  flow 
channel  because  the  relative  humidity  condition  in  the  gas  flow 
channels  is  the  same.  Water  vapor  mass  fraction  shows  a  tendency 
to  increase  as  the  current  density  increases.  This  can  be  explained 
by  looking  at  the  oxygen  contours  of  Fig.  11(b).  Oxygen  concen¬ 
tration  decreases  as  current  density  increases  due  to  the  increasing 
consumption  rate,  thus,  due  to  the  non-reactivity  of  nitrogen  then 
the  water  vapor  concentration  must  increase. 

The  parametric  studies  presented  so  far  on  saturation  and 
species  distribution  have  shown  some  important  capabilities  of  the 


Fig.  12.  Parametric  study  of  water  vapor  distribution  in  the  GDL  (a)  Cell  operating  at  1.4  A  cm  2.  (b)  Air  inlet  relative  humidity  of  95%. 
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Fig.  13.  Parametric  study  of  the  net  phase  change  rate  in  the  GDL.  (a)  Cell  operating  at  1.4  A  cm  2.  (b)  Air  inlet  relative  humidity  of  95%. 


model  developed.  However,  results  like  the  ones  presented  in 
Fig.  11  would  be  better  explained  by  looking  at  the  phase  change 
rates  in  the  computational  domain.  Furthermore,  it  is  interesting  to 
look  at  how  the  phases  interact  inside  the  GDL.  Fig.  13  shows 
a  parametric  study  of  the  phase  change  rate  in  the  GDL.  According 
to  the  phase  change  model  (Eq.  (10)),  the  phase  change  driving 
mechanism  is  the  non-equilibrium  between  local  vapor  pressure 
and  the  saturated  vapor  pressure  at  the  local  temperature.  The 
convention  of  the  present  modeling  efforts  is  that  a  positive  phase 
change  rate  represents  a  gain  in  water  vapor  (evaporation)  and 
a  negative  phase  change  rate  represents  condensation.  It  can  be 
observed  in  Fig.  13(a)  that  the  higher  evaporation  rate  is  for 
a  condition  of  60%  relative  humidity  and  the  tendency  is  to  reduce 
the  evaporation  rate  as  the  inlet  relative  humidity  condition 
increases.  The  curves  depicted  in  Fig.  13  show  a  slight  tendency  to 
reduce  the  magnitude  of  evaporation  rate  in  the  regions  close  to  the 
catalyst  layer  interfaces;  such  a  behavior  is  due  to  the  increase  in 
relative  humidity  near  the  catalyst  layer  due  to  the  higher  content 
of  liquid  water.  Another  important  feature  observed  in  the  curves  of 
Fig.  13  is  that  they  resemble  the  saturation  profiles  shown  in  Fig.  10. 


This  is  because  the  phase  change  model  depends  upon  the  exis¬ 
tence  of  saturation  in  order  to  have  evaporation.  Condensation  can 
occur  at  regions  saturated  with  liquid  water  but  only  if  the  relative 
humidity  is  higher  than  1. 

The  results  shown  in  Fig.  13(a)  predict  that  only  evaporation 
occurs  at  the  operating  conditions  studied.  Fig.  13(b)  shows  how 
the  model  evolves  from  evaporation  to  condensation  at  a  fixed  high 
inlet  humidity  condition  (95%)  while  current  density  is  increased. 
At  the  lower  current  density  condition,  0.16  A  cm-2,  the  phenom¬ 
enon  inside  the  GDL  is  purely  evaporative;  however,  as  the  current 
density  condition  increases,  liquid  water  saturation  and  local 
relative  humidity  increases  from  the  gas  flow  channels  interface  to 
the  catalyst  layer  interface.  These  phenomena  create  the  peak 
shown  in  the  curves  and  then  the  decrease  of  phase  change  rate 
to  a  condensation  condition  at  the  highest  current  density 
condition. 

Local  relative  humidity  profiles  were  obtained  from  the  para¬ 
metric  studies  shown  previously.  Fig.  14  shows  profiles  of  local 
relative  humidity  at  different  operating  conditions.  Fig.  14(a) 
shows  a  tendency  to  increase  the  local  relative  humidity  from  the 


Fig.  14.  Parametric  study  of  local  relative  humidity  in  the  GDL  (a)  Cell  operating  at  1.4  A  cm  2.  (b)  Air  inlet  relative  humidity  of  95%. 
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gas  flow  channel  interface  to  the  catalyst  layer  interface.  Such  an 
increase  in  relative  humidity  is  most  likely  attributed  to  the 
consumption  of  oxygen  and  the  consequent  increase  of  concen¬ 
tration  of  water  vapor.  The  increase  of  relative  humidity  is  also 
related  to  evaporation  (see  Fig.  13(a))  due  to  the  high  driving 
mechanism  for  this  phase  change  phenomenon  at  low  humidity 
conditions  in  the  gas  flow  channels.  The  amount  of  relative 
humidity  gain  is  reduced  as  the  gas  flow  channel  humidity 
condition  increases.  Fig.  14(b)  shows  that  local  relative  humidity 
increases  in  the  vicinity  of  the  catalyst  layers  as  the  average  current 
density  increases.  This  is  explained  because  oxygen  is  depleted 
faster  at  higher  current  densities,  water  vapor  concentration 
increases,  and  therefore  relative  humidity  increases  too.  Due  to  the 
high  humidity  condition  in  the  gas  flow  channels,  the  phase 
change  rate  that  depends  upon  relative  humidity  gradients  turns 


a  A  0.167 


from  evaporation  to  condensation,  see  Fig.  13(b).  That  phenom¬ 
enon  can  be  explained  by  looking  at  Fig.  14(b).  It  is  observed  how 
at  2.1  A  cm-2,  the  local  relative  humidity  is  greater  than  1,  hence, 
the  condensation  phenomenon  substitutes  evaporation  near  the 
catalyst  layer  interface. 


4.2.  2-D  model  results  (validation ) 

It  is  commonly  observed  that  the  objective  of  most  of  the 
numerical  models  developed  for  fuel  cells  is  to  predict  polarization 
behavior.  Currently,  models  reported  in  the  literature  generally  use 
a  two-phase  approach  and  most  of  them,  report  to  predict  polari¬ 
zation  adequately  although  the  models  are  based  on  assumed 
properties  of  the  GDL  material  and  also  assumed  electrochemical 
parameters.  The  novelty  of  the  model  presented  in  this  work  is  to 
predict  polarization  by  using  actual  experimentally  determined 
properties  of  the  GDL  material  and  experimentally  obtained  elec¬ 
trochemical  properties  as  well.  In  Section  3.2  two  polarization 
curves  were  obtained  for  a  cell  assembled  using  non-wet-proofed 
and  wet-proofed  GDL  material.  The  same  MEA  material  was 
taken  from  the  same  stock  and  the  cells  were  tested  at  the  same 
operating  conditions. 

The  main  objective  of  this  work  is  to  determine  if  polarization 
curves  can  be  predicted  by  using  experimentally  determined  water 
transport  properties  of  the  GDL  material,  given  that  many  of  the 
most  significant  losses  occur  within  this  region  of  the  cell,  namely, 
diffusion  of  species,  liquid  water  saturation  of  electrodes,  and 
electric  charge  transfer.  This  hypothesis  was  tested  and  the  steps 
followed  are  summarized  here:  (1)  two  cells  were  assembled  using 
different  GDL  material  (treated  and  untreated  with  PTFE)  and  in 
such  a  way  that  a  compression  of  ~20%  was  obtained  for  the  GDL 
material;  (2)  the  water  transport  characteristics  of  GDL  materials  at 
compressions  of  ~20%  were  obtained  experimentally  [36];  (3)  data 
from  the  GDLs  and  the  Tafel  plot  (electrochemical  parameters) 
were  implemented  in  a  numerical  model  of  the  air-cathode  of 
a  PMEFC  to  obtain  the  polarization  curves  of  each  cell.  Fig.  15 
presents  a  comparison  between  experimental  polarization  curves 
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Fig.  16.  Saturation  contours  of  liquid  water  saturation,  (a)  Non-wet-proofed  GDL  (b)  Wet-proofed  GDL 
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and  the  polarization  curves  predicted  using  the  numerical  2-D 
model  developed.  It  can  be  clearly  seen  that  a  good  prediction  of 
polarization  was  obtained  using  the  numerical  model.  The  three 
overpotential  zones  are  observed  in  the  numerical  model  results; 
however,  the  concentration  zone  is  where  both  numerical  results, 
for  the  untreated  and  the  wet-proofed  GDL  model,  differ  with  the 
experimental  data.  Actually,  this  region  was  very  difficult  to 
calculate  beyond  an  operating  condition  of  0.4  V  due  to  divergence 
of  the  solution  process. 

Under  high  current  density  conditions  (e.g.,  >  0.9  A  cm-2  )  the 
concentration  loss  effects  become  more  significant  for  both  fuel 
cells.  At  this  high  current  density  condition  liquid  water  is  produced 
at  high  rates  and  has  partially  saturated  the  electrode.  As  seen  in 
the  following  section,  saturation  is  particularly  high  under  the 
shoulders  of  the  collector  plates.  Sole  [42]  has  suggested  that  this 
area  under  the  collector  plate  shoulder  exhibits  the  lowest  contact 
resistance  and  that  excessive  saturation  in  a  region  that  would 
otherwise  produce  the  highest  current,  has  a  more  detrimental 
effect  than  is  predicted  by  models  which  assume  a  uniform  contact 
resistance  over  the  surface  of  the  MEA.  This  rationale  may  explain 
why  models  such  as  the  one  presented  here  (which  does  not 
address  variations  in  contact  resistance)  overpredict  cell  perfor¬ 
mance  at  very  high  current  density. 

Fig.  16  presents  liquid  water  saturation  contours  for  both  cells, 
using  untreated  and  treated  GDL  materials  at  an  operating  condi¬ 
tion  of  1  A  cm-2.  The  saturation  contours  for  the  different  materials 
are  shown  using  different  scale  due  to  the  considerable  differences 
between  both  saturation  levels.  It  can  be  observed  that  the  cell 
using  non-wet-proofed  material  presents  more  saturation  than  the 
cell  using  wet-proofed  material,  such  a  result  was  expected  and 
makes  sense.  In  spite  of  the  differences  between  the  magnitudes  of 
saturation  levels,  similar  characteristics  can  be  observed  and  will  be 
discussed  here.  It  is  observed  in  both  cases  shown  in  Fig.  16  that  the 
lowest  saturation  levels  occur  near  the  interface  with  the  gas  flow 
channels.  This  is  evident  because  at  that  interface,  the  gas  stream 
removes  liquid  water  by  a  convective  mechanism.  It  is  also 
observed  that  the  highest  liquid  water  saturation  levels  are  present 
in  the  region  below  the  land  (under  the  current  collectors).  Such  an 
accumulation  of  water  can  be  explained  due  to  the  impermeability 
of  the  current  collector  to  the  water,  and  due  to  the  long  pathway 
between  the  water  generated  below  the  land  and  the  interface  with 
the  gas  flow  channels. 

5.  Conclusions 

Much  of  the  continuum  modeling  approaches  for  PEM  fuel  cells 
reported  in  the  literature  use  empirical  correlations  for  materials 
like  well-sorted  sand  as  inputs  to  the  water  transport  equation. 
Only  recently  has  research  been  conducted  and  reported  in  the 
literature  in  regards  to  the  characterization  of  GDL  material  in  order 
to  provide  actual  water  transport  properties  such  as  capillary 
pressure,  porosity,  and  permeability.  In  the  present  work  such 
properties  were  taken  from  experimentally  reported  data  for 
carbon  paper  Toray  090  [36]. 

Experimental  performance  was  obtained  for  fuel  cells  using 
untreated  GDL  material  and  also  using  treated  GDL  material  with 
20%  of  PTFE.  It  is  well  known  that  wet-proofing  is  applied  to  increase 
performance  of  the  cell  by  avoiding  electrode  flooding.  However,  it 
is  also  reported  in  the  literature  that  the  amount  of  wet-proofing 
must  be  optimized  for  the  intended  operating  conditions.  In  the 
experiments  conducted  here  (at  low  humidity  conditions)  it  was 
observed  that  the  cell  using  an  untreated  GDL  performed  better 
than  the  cell  using  a  wet-proofed  GDL.  These  results  were  attributed 
to  the  (1 )  reduction  of  the  pore  space  due  to  the  PTFE  content  in  the 
GDL;  (2)  creation  of  a  more  tortuous  path  due  to  the  PTFE  treatment; 


(3)  reduction  of  the  GDL  conductivity  due  to  the  PTFE  sintered  in  the 
carbon  fibers  of  the  GDL;  and  (4)  reduction  of  MEA  humidification 
due  to  the  hydrophobic  GDL  used  in  the  anode  and  cathode  side. 
From  these  experiments,  electrochemical  parameters  such  as  the 
transfer  coefficient,  the  exchange  current  density,  and  the  operating 
cell  voltage  were  determined  for  the  MEA  material  used. 

Experimental  characterization  data  of  GDL  material  and  elec¬ 
trochemical  parameters  determined  for  the  MEA  material  were 
used  as  inputs  of  two  numerical  models  of  a  PEMFC.  A  1-D  model 
was  developed  in  order  to  assess  the  effect  of  using  experimental 
water  transport  characteristics  of  a  GDL  on  liquid  water  saturation 
distribution.  The  liquid  water  saturation  profiles  obtained  using 
experimentally  determined  properties  of  the  GDL  were  compared 
to  those  calculated  using  the  most  commonly  used  empirical 
correlations  reported  in  the  literature.  It  was  found  that  the 
empirical  correlations  underpredict  saturation  levels  in  the  GDL. 
The  effect  of  wet-proofing  was  also  studied  and  a  remarkable 
reduction  of  saturation  level  was  observed  at  the  operating 
conditions  simulated.  A  parametric  study  was  conducted  to 
demonstrate  the  reliability  of  the  numerical  model  to  predict 
species  distribution  and  water  saturation.  The  results  obtained 
from  this  parametric  study  confirmed  an  appropriate  description  of 
the  physics  governing  the  operation  of  PEM  fuel  cells. 

At  the  end  and  as  a  closing  study,  a  2-D  model  was  developed 
seeking  to  replicate  the  polarization  curves  obtained  using  different 
GDL  material.  It  was  observed  that  by  varying  the  water  transport 
characteristics  of  the  GDL,  a  good  agreement  between  experimental 
and  numerical  data  was  attained  for  current  density  as  high  as  1  A 
cm-2.  With  these  results,  it  was  demonstrated  that  polarization  of 
PEM  fuel  cells  using  different  GDL  material  can  be  predicted  just  by 
using  adequate  GDL  and  water  transport  properties.  Further 
investigation  of  the  capabilities  of  the  numerical  model  was  con¬ 
ducted  by  obtaining  saturation  distribution  in  the  regions  below 
the  current  collectors  and  under  the  gas  flow  channels.  It  was  found 
that  the  liquid  water  saturation  is  higher  under  the  land  of  the 
current  collectors  than  under  the  channels.  Such  results  were  ex¬ 
pected  and  match  with  other  works  reported  in  the  literature. 
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